
# Replication file for Figure 1 (Body); Appendix G.1.: Tables 11-16

#install.packages("Hmisc")
require(Hmisc)
#install.packages("stargazer")
require(stargazer)

load("Brexit_Rep.RData")

leavers.sub <- brexit.sub2[which(brexit.sub2$leaver == 1),]
remainers.sub <- brexit.sub2[which(brexit.sub2$leaver == 0),]

dvs <- Cs(migs.take.jobs.6,migs.more.terror.6,closed.immigrants.6,hurt.standing.refugee.6,
          threaten.culture.refugee.6,overwhelm.welfare.refugee.6)

#### Controls (imputed income): Tables 11 and 12 ####

right.side.eq <- "~ post + Gender + Age + Hincome.imp + socialgrade + Education + employed + children + marstat + as.factor(GOR_full)"

c.ate.l <- vector(length(dvs), mode = "list")
names(c.ate.l) <- dvs

c.ate.r <- vector(length(dvs), mode = "list")
names(c.ate.r) <- dvs


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  c.ate.l[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = leavers.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
  print(summary(c.ate.l[[dvs[i]]]))
  
}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  c.ate.r[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = remainers.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
  print(summary(c.ate.r[[dvs[i]]]))
  
}


varlabels <- c("Treat", "Male", "Age", "Income (Imp.)", "Social Grade", "Education","Employed", "Children (=1)", "Married", "North West", 
               "Yorkshire and the Humber", "East Midlands", "West Midlands", "East of England","London", "South East", 
               "South West", "Wales", "Scotland", "Constant")

#Table 11
stargazer(c.ate.l[[dvs[1]]], c.ate.r[[dvs[1]]],c.ate.l[[dvs[2]]],c.ate.r[[dvs[2]]],c.ate.l[[dvs[3]]],c.ate.r[[dvs[3]]], covariate.labels = varlabels)
#Table 12
stargazer(c.ate.l[[dvs[5]]], c.ate.r[[dvs[5]]],c.ate.l[[dvs[6]]],c.ate.r[[dvs[6]]],c.ate.l[[dvs[4]]],c.ate.r[[dvs[4]]], covariate.labels = varlabels)



### Controls (non-imputed income): #### Tables 13 and 14  #####

cni.ate.l <- vector(length(dvs), mode = "list")
names(cni.ate.l) <- dvs

cni.ate.r <- vector(length(dvs), mode = "list")
names(cni.ate.r) <- dvs

right.side.eq <- "~ post + Gender + Age + Hincome.rec + socialgrade + Education + employed + children + marstat + as.factor(GOR_full)"


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  cni.ate.l[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = leavers.sub, 
                                          weights = w8w6), list(.modelformula = modelformula)))
  
  print(summary(cni.ate.l[[dvs[i]]]))
  
}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  cni.ate.r[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = remainers.sub, 
                                          weights = w8w6), list(.modelformula = modelformula)))
  
  print(cni.ate.r[[dvs[i]]])
  
}

varlabels <- c("Treat", "Male", "Age", "Income (Orig.)", "Social Grade", "Education", "Employed", "Children (=1)", "Married", "North West", 
               "Yorkshire and the Humber", "East Midlands", "West Midlands", "East of England","London", "South East", 
               "South West", "Wales", "Scotland", "Constant")


#Table 13
stargazer(cni.ate.l[[dvs[1]]], cni.ate.r[[dvs[1]]],cni.ate.l[[dvs[2]]],cni.ate.r[[dvs[2]]],cni.ate.l[[dvs[3]]],cni.ate.r[[dvs[3]]], covariate.labels = varlabels)

#Table 14
stargazer(cni.ate.l[[dvs[5]]], cni.ate.r[[dvs[5]]],cni.ate.l[[dvs[6]]],cni.ate.r[[dvs[6]]],cni.ate.l[[dvs[4]]],cni.ate.r[[dvs[4]]], covariate.labels = varlabels)


### Controls (minus income) #### Tables 15 and 16

cmi.ate.l <- vector(length(dvs), mode = "list")
names(cmi.ate.l) <- dvs

cmi.ate.r <- vector(length(dvs), mode = "list")
names(cmi.ate.r) <- dvs

right.side.eq <- "~ post + Gender + Age + socialgrade + Education + employed + children + marstat + as.factor(GOR_full)"


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  cmi.ate.l[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = leavers.sub, 
                                          weights = w8w6), list(.modelformula = modelformula)))
  
  print(summary(cmi.ate.l[[dvs[i]]]))
  
}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  cmi.ate.r[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = remainers.sub, 
                                          weights = w8w6), list(.modelformula = modelformula)))
  
  print(summary(cmi.ate.r[[dvs[i]]]))
  
}

varlabels <- c("Treat", "Male", "Age", "Social Grade", "Education","Employed", "Children (=1)", "Married", "North West", 
               "Yorkshire and the Humber", "East Midlands", "West Midlands", "East of England","London", "South East", 
               "South West", "Wales", "Scotland", "Constant")


#### Tables 15 and 16

stargazer(cmi.ate.l[[dvs[1]]], cmi.ate.r[[dvs[1]]],cmi.ate.l[[dvs[2]]],cmi.ate.r[[dvs[2]]],cmi.ate.l[[dvs[3]]],cmi.ate.r[[dvs[3]]], covariate.labels = varlabels)

stargazer(cmi.ate.l[[dvs[5]]], cmi.ate.r[[dvs[5]]],cmi.ate.l[[dvs[6]]],cmi.ate.r[[dvs[6]]],cmi.ate.l[[dvs[4]]],cmi.ate.r[[dvs[4]]], covariate.labels = varlabels)



###### Figure 1 (body) ##########

cis.l <- matrix(NA,length(dvs),2)
b.l <- matrix(NA,length(dvs),1)
cis.r <- matrix(NA,length(dvs),2)
b.r <- matrix(NA,length(dvs),1)

for (i in 1:length(dvs))
{
  b.l[i] <- c.ate.l[[i]]$coefficients[2]
  print(c.ate.l[[i]]$coefficients[2])
  ci <- confint(c.ate.l[[i]])
  #print(ci)
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.l[i,] <- cbind(lb,ub)
}

for (i in 1:length(dvs))
{
  b.r[i] <- c.ate.r[[i]]$coefficients[2]
  ci <- confint(c.ate.r[[i]])
  print(b.r[i])
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.r[i,] <- cbind(lb,ub)
}


y.lab <- seq(from = 1, to=length(dvs),by=1)
vnames <-c(migs.take.jobs.6 = "Migs. Take\nJobs",migs.more.terror.6 = "Migs. Bring\nTerror", closed.immigrants.6 = "Not\nOpen to\nMigs.",
           hurt.standing.refugee.6 = "Refs. Don't\nImprove\nUK Image", threaten.culture.refugee.6 = "Refs.\nThreaten\nCulture", overwhelm.welfare.refugee.6 ="Refs.\nOverwhelm\nServices")


jitter <-  0.1

pdf("Fig1.pdf", height = 5, width=5)
par(mar =c(2,1.1, 1.1, 0.01), oma=c(1,5,.2,.2), cex.axis=.8, cex.main=.9)
plot(b.l, y.lab, type="n", ylab ="", xlab = "", yaxt="n", xlim = c(-0.5,0))
axis(2, at=y.lab, labels=vnames, las = 2, cex.axis=0.8)
mtext("Difference in Means",side=1,line=2,outer=F, cex = 0.8)
#segments(br.take.jobs.r$coefficients[2],x.lab, m.ci.2.lp, yr.lab, lty=1)
points(b.l, y.lab, pch= 19, cex=.8, col = "black")
segments(cis.l[,1], y.lab, cis.l[,2], y.lab, lty=1, col = "black")
segments(cis.r[,1], y.lab+jitter, cis.r[,2], y.lab+jitter, lty=1, col = "black")
points(b.r, y.lab+jitter, pch= 19, cex=.8, col = "white")
points(b.r, y.lab+jitter, pch= 1, cex=.8)
abline(v=0)
legend(-0.5, 6, pch=c(19,1), legend=c("Leave", "Remain"), cex =.8, border = "black")
dev.off()

